clc;  
%——————————低通滤波器滤除肌电信号——————————  



%采样获得的数据
data=importdata('data4.txt');
M=data';

TIME=linspace(0,18.67,1868);%时间
Fs=100;                        %采样频率  
fp=10;
fs=15;                    %通带截止频率，阻带截止频率  

rp=1.4;
rs=1.6;                    %通带、阻带衰减  
wp=2*pi*fp;
ws=2*pi*fs;     

[n,wn]=buttord(wp,ws,rp,rs,'s');     %’s’是确定巴特沃斯模拟滤波器阶次和3dB截止模拟频率  
[z,P,k]=buttap(n);   %设计归一化巴特沃斯模拟低通滤波器，z为极点，p为零点和k为增益  
[bp,ap]=zp2tf(z,P,k)  %转换为Ha(p),bp为分子系数，ap为分母系数  
[bs,as]=lp2lp(bp,ap,wp) %Ha(p)转换为低通Ha(s)并去归一化，bs为分子系数，as为分母系数  
  
[hs,ws]=freqs(bs,as);         %模拟滤波器的幅频响应  
[bz,az]=bilinear(bs,as,Fs);     %对模拟滤波器双线性变换  
[h1,w1]=freqz(bz,az);         %数字滤波器的幅频响应  
m=filter(bz,az,M);  
  
figure  
freqz(bz,az);
title('巴特沃斯低通滤波器幅频曲线');  
        
figure  
subplot(2,1,1);  
%plot(TIME,M(:,1));  
plot(TIME,M);  
xlabel('t(s)');
ylabel('mv');
title('原始心电信号波形');
grid;  
  
subplot(2,1,2);  
plot(TIME,m);  
xlabel('t(s)');ylabel('mv');
title('低通滤波后的时域图形');
grid;  
     
% -----------------------频谱图-----------------------------------

N=512  %这个表示的是512个傅里叶谐波
n=0:N-1;  
mf=fft(M,N);               %进行频谱变换（傅里叶变换）  
mag=abs(mf);  
f=(0:length(mf)-1)*Fs/length(mf);  %进行频率变换  
  
figure  
subplot(2,1,1)  
plot(f,mag);axis([0,150,1,1500]);grid;      %画出频谱图  
xlabel('频率(HZ)');
ylabel('幅值');
title('心电信号频谱图'); 
grid;
  
mfa=fft(m,N);                    %进行频谱变换（傅里叶变换）  
maga=abs(mfa);  
fa=(0:length(mfa)-1)*Fs/length(mfa);  %进行频率变换  
subplot(2,1,2)  
plot(fa,maga);axis([0,150,1,1500]);grid;  %画出频谱图  
xlabel('频率(HZ)');
ylabel('幅值');
title('低通滤波后心电信号频谱图');  
grid;
   % -----------------------心电信号的功率谱-----------------------------------   
wn=M;  
P=10*log10(abs(fft(wn).^2)/N);  
f=(0:length(P)-1)/length(P);  
figure  
plot(f,P);
grid  
xlabel('归一化频率');
ylabel('功率(dB)');
title('心电信号的功率谱');  




%--------------------------检查下两个变量维度------------------------------------
m2=m;
size(m2)
size(TIME)


%------------------IIR零相移数字滤波器纠正基线漂移-------------------
% Wp=1.4*2/Fs;     %通带截止频率 
% Ws=0.6*2/Fs;     %阻带截止频率 

Wp=1.4*2/Fs;     %通带截止频率 
Ws=0.6*2/Fs;     %阻带截止频率 
devel=0.005;    %通带纹波 
Rp=20*log10((1+devel)/(1-devel));   %通带纹波系数  
Rs=20;                          %阻带衰减 
[N Wn]=ellipord(Wp,Ws,Rp,Rs,'s');   %求椭圆滤波器的阶次 
[b a]=ellip(N,Rp,Rs,Wn,'high');       %求椭圆滤波器的系数 
[hw,w]=freqz(b,a,512);   
result =filter(b,a,m2); 
 
figure
freqz(b,a);
figure
subplot(211); 
plot(TIME,m2); 
xlabel('t(s)');
ylabel('幅值');
title('Butterworth低通滤波器输出信号');
axis([0,20,-100,500]);
grid




subplot(212); 
plot(TIME,result); 
xlabel('t(s)');
ylabel('幅值');
title('基线噪声零相移滤波后信号');
axis([0,20,-100,150]);
grid
  
figure
N=512;

subplot(2,1,1);plot(abs(fft(m2))*2/N);
xlabel('频率(Hz)');
ylabel('幅值');
title('Butterworth低通滤波器输出信号频谱');
axis([0,2000,1,40]);
grid;



subplot(2,1,2);
plot(abs(fft(result))*2/N);
xlabel('频率(Hz)');
ylabel('幅值');
title('基线噪声零相移滤波后');
grid;

subplot(2,1,2);
plot(abs(fft(result))*2/N);
xlabel('基线噪声零相移滤波后信号频谱');
ylabel('幅值');
grid;





%----------------------------------------------------------------------------------------
sprintf('-----------findpeaks----------------------')
plot(TIME,result)
grid;
hold on
% [pks,locs] = findpeaks(result)

Threshold = (max(result) - min(result))*0.4 + min(result);
sprintf('-----------Threshold---------------------')
max(result)
min(result)
Threshold

% result
[pks,locs] =findpeaks(result,100,'MinPeakDistance',0.8,'MinPeakHeight',Threshold);
% 这里因为已经填写了100Hz的采样频率，所以findpeaks中的全程Distance是1868*1/100=18.68
% 所以这里的18.68成为了MinPeakDistance的范围上限

sprintf('-----------pks locs---------------------')
pks
locs
% plot(locs*0.01,pks);

sprintf('-----------下面计算心率---------------------')
locs_valid = [locs(3:end)];%去掉开头两条异常的，
delta_time=diff(locs_valid)%每两次心跳时间间隔序列，剩下16条，计算间隔差，所以是15条数据
heart_rate=60/mean(delta_time)%60/心跳一次平均消耗的时间





sprintf('-----------经过两次滤波后的数据保存到硬盘---------------------')
 
fid=fopen('result.txt','w');
% fprintf(fid,'%7s %9s\r\n','Radius','Area');%写入列名
fprintf(fid,'%8d\r\n',result);%写入每列的数据
fclose(fid);


sprintf('-----------经过两次滤波后的数据保存到硬盘---------------------')

fid=fopen('locs.txt','w');
% fprintf(fid,'%7s %9s\r\n','Radius','Area');%写入列名
fprintf(fid,'%8d\r\n',locs);%写入每列的数据
fclose(fid);


